#functino to compute distance from targeted moments
function Compute_error(prim_2000::Primitives, output_simul::Vector{DenseArray{Float64, N} where N}, package_moments::Array{Array{Float64,2},1}; report::Int64 = 1)
    @unpack nl, na, divs, prox_mat, pops = prim_2000 #unpack
    data_simul, moving_mat_simul = output_simul[1], output_simul[2] #unpack simualted output
    moments, moving_mat = package_moments[1], package_moments[2] #unpack target moments
    nrows, ncol = length(data_simul[:,1]), length(data_simul[1,:]) #useful for later
    attendance_coll, attendance_noncoll = package_moments[3], package_moments[4]
    moments_race = package_moments[7]
    #error = 0 #preallocate

    err_vec = zeros(234) #number of moments

    #####separete into college, non-college; construct weights
    data_simul_hs = data_simul[data_simul[:,13].==1, :]
    data_simul_coll = data_simul[data_simul[:,13].==2, :]
    weights_all = fweights(data_simul[:,14]) #define population weights
    weights_hs = fweights(data_simul_hs[:,14])
    weights_coll = fweights(data_simul_coll[:,14])

    ###########IGE###########
    ige = cor(data_simul[:,ncol-1], data_simul[:,ncol])
    err_vec[1] = ((moments[1] - ige)/moments[1]) #add to error term

    ###########Earnings Moments###########

    ####High School
    mean_earn_hs_p2 = mean(data_simul_hs[:,8], weights(weights_hs))
    std_earn_hs_p2 = std(data_simul_hs[:,8], weights(weights_hs); corrected = false)
    mean_earn_hs_growth = mean(log.(data_simul_hs[:,9]) - log.(data_simul_hs[:,8]), weights(weights_hs))
    std_earn_hs_growth = std(log.(data_simul_hs[:,9]) - log.(data_simul_hs[:,8]), weights(weights_hs); corrected=false)
    mean_earn_hs_p3 = mean(data_simul_hs[:,9], weights(weights_hs))
    std_earn_hs_p3 = std(data_simul_hs[:,9], weights(weights_hs); corrected=false)

    err_vec[2] = ((moments[2] - mean_earn_hs_p2)/moments[2])
    err_vec[3] = ((moments[3] - std_earn_hs_p2)/moments[3])
    err_vec[4] = ((moments[6] - mean_earn_hs_growth)/moments[6])
    err_vec[5] = ((moments[7] - std_earn_hs_growth)/moments[7])
    err_vec[6] = ((moments[10] - mean_earn_hs_p3)/moments[10])
    err_vec[7] = ((moments[11] - std_earn_hs_p3)/moments[11])

    ####College
    mean_earn_coll_p2 = mean(data_simul_coll[:,8], weights(weights_coll))
    std_earn_coll_p2 = std(data_simul_coll[:,8], weights(weights_coll); corrected=false)
    mean_earn_coll_growth = mean(log.(data_simul_coll[:,9]) - log.(data_simul_coll[:,8]), weights(weights_coll))
    std_earn_coll_growth = std(log.(data_simul_coll[:,9]) - log.(data_simul_coll[:,8]), weights(weights_coll); corrected=false)
    mean_earn_coll_p3 = mean(data_simul_coll[:,9], weights(weights_coll))
    std_earn_coll_p3 = std(data_simul_coll[:,9], weights(weights_coll); corrected=false)

    err_vec[8] = ((moments[4] - mean_earn_coll_p2)/moments[4])
    err_vec[9] = ((moments[5] - std_earn_coll_p2)/moments[5])
    err_vec[10] = ((moments[8] - mean_earn_coll_growth)/moments[8])
    err_vec[11] = ((moments[9] - std_earn_coll_growth)/moments[9])
    err_vec[12] = ((moments[12] - mean_earn_coll_p3)/moments[12])
    err_vec[13] = ((moments[13] - std_earn_coll_p3)/moments[13])

    #println("Error so far: ", error)

    #######average time investment#######
    time_inv_avg = mean(data_simul[:,7], weights(weights_all))
    err_vec[14] = ((moments[23] - time_inv_avg)/moments[23])
    index = 15
    coll_weight = 0.34

    ##########COLLEGE ATTENDANCE
    data_simul[:,13].-=1 #change college attainment to boolean; makes computing means easier

    ##parents without college degree
    data_simul_nonpcoll = data_simul[data_simul[:,12].==1,:] #parents who did not have a college degree
    coll_p_noncoll = mean(data_simul_nonpcoll[:,13])
    attendance_noncoll_simul = zeros(4, na)

    for i = 1:4, j = 1:na

        #restrict data to relevant ability level
        restrict_ability = data_simul_nonpcoll[data_simul_nonpcoll[:,11].==j, :]

        #restrict to parent income quartile
        minpctile = (i-1)*25
        maxpctile = (i)*25
        restrict_pctile_pt1 = restrict_ability[restrict_ability[:,ncol-1].>=minpctile, :]
        restrict_pctile_pt2 = restrict_pctile_pt1[restrict_pctile_pt1[:,ncol-1].<maxpctile, :]


        if length(restrict_pctile_pt2[:,1]) ==0
            println("awooga", i, j)
        end

        weights_temp = fweights(restrict_pctile_pt2[:,14])
        attendance_noncoll_simul[i, j] = mean(restrict_pctile_pt2[:,13], weights(weights_temp))
    end

    ##parents with college degree
    data_simul_pcoll = data_simul[data_simul[:,12].==2,:] #parents who did not have a college degree
    coll_p_coll = mean(data_simul_pcoll[:,13])
    attendance_coll_simul = zeros(4, na)

    for i = 1:4, j = 1:na

        #restrict data to relevant ability level
        restrict_ability = data_simul_pcoll[data_simul_pcoll[:,11].==j, :]

        #restrict to parent income quartile
        minpctile = (i-1)*25
        maxpctile = (i)*25
        restrict_pctile_pt1 = restrict_ability[restrict_ability[:,ncol-1].>=minpctile, :]
        restrict_pctile_pt2 = restrict_pctile_pt1[restrict_pctile_pt1[:,ncol-1].<maxpctile, :]

        if length(restrict_pctile_pt2[:,1]) ==0
            println("awooga2", i, j)
        end

        weights_temp = fweights(restrict_pctile_pt2[:,14])
        attendance_coll_simul[i, j] = mean(restrict_pctile_pt2[:,13], weights(weights_temp))
    end

    for i = 1:4, j = 1:na
        err_vec[index] = ((attendance_noncoll[i, j] - attendance_noncoll_simul[i, j])/coll_weight)/sqrt(40)
        index += 1
    end

    for i = 1:4, j = 1:na
        if attendance_coll[i, j] != 999 #flag for being fewer than 25 observations
            err_vec[index]  = ((attendance_coll[i, j] - attendance_coll_simul[i, j])/coll_weight)/sqrt(40)
        end
        index += 1
    end

    ###########Migration Moments#############

    #W = BlockDiagonal([mat_mig_moments, mat_mig_dynamic, mat_mig_mat])

    ####Migration Rates for College Grads, Non-Grads
    mig_hs = mean(data_simul_hs[:,15], weights(weights_hs))
    mig_coll = mean(data_simul_coll[:,15], weights(weights_coll))
    err_vec[55] = ((moments[14] - mig_hs)/moments[14])
    err_vec[56] = ((moments[15] - mig_coll)/moments[15])

    ####HC Differences Between Grads, Non-Grads

    #hs
    restrict_move_hs = data_simul_hs[data_simul_hs[:,1].!=data_simul_hs[:,2],:]
    restrict_stay_hs = data_simul_hs[data_simul_hs[:,1].==data_simul_hs[:,2],:]
    hc_index = 10
    weights_move_hs = fweights(restrict_move_hs[:,14])
    weights_stay_hs = fweights(restrict_stay_hs[:,14])
    hc_move_hs = mean(restrict_move_hs[:,hc_index], weights(weights_move_hs))
    hc_stay_hs = mean(restrict_stay_hs[:,hc_index], weights(weights_stay_hs))
    err_vec[57] = ((moments[16] - (hc_move_hs-hc_stay_hs))/moments[16])

    #coll
    restrict_move_coll = data_simul_coll[data_simul_coll[:,1].!=data_simul_coll[:,2],:]
    restrict_stay_coll = data_simul_coll[data_simul_coll[:,1].==data_simul_coll[:,2],:]
    hc_index = 10
    weights_move_coll = fweights(restrict_move_coll[:,14])
    weights_stay_coll = fweights(restrict_stay_coll[:,14])
    hc_move_coll = mean(restrict_move_coll[:,hc_index], weights(weights_move_coll))
    hc_stay_coll = mean(restrict_stay_coll[:,hc_index], weights(weights_stay_coll))
    err_vec[58] = ((moments[17] - (hc_move_coll-hc_stay_coll))/moments[17])

    ###average
    mig_avg = mean(diag(moving_mat_simul))
    #err_vec[59] = ((moments[18]-mig_avg)/moments[18])

    ###stdev
    mig_stdev = std(diag(moving_mat_simul); corrected=false)
    #err_vec[60] = ((moments[19] - mig_stdev)/moments[19])

    ###correlation between population and absorption
    state_receipts = zeros(50)
    for l = 1:50
        state_receipts[l] = sum(moving_mat_simul[:,l]) - diag(moving_mat_simul)[l] #summing migration probabiltiies, leaving stayers out
    end
    corr_pop_mig_in = cor(state_receipts, pops)
    #error += ((moments[20] - corr_pop_mig_in)/moments[20])^2

    ####correlation between population and staying
    corr_pop_stay = cor(pops, diag(moving_mat_simul))
    #error += ((moments[21] - corr_pop_stay)/moments[21])^2


    ###share of moves to proximal state
    prox_move_mat = moving_mat_simul.*prox_mat
    prox_move_rate = (sum(prox_move_mat) - sum(diag(prox_move_mat)))/nl
    mig_rate = (sum(moving_mat_simul) - sum(diag(moving_mat_simul)))/nl
    prox_share = prox_move_rate/mig_rate
    err_vec[61] = ((moments[22]-prox_share)/moments[22])

    ######dynamic migration stuff
    mig_p2_rate = mean(data_simul[:,17], weights(weights_all))
    #mig_home_rate = mean(data_simul[:,18], weights(weights_all))

    err_vec[62] = ((moments[26] - mig_p2_rate)/moments[26])
    #error += ((moments[27] - mig_home_rate)/moments[27])^2


    ##overall correlation outflow
    #=
    mig_corr = cor(diag(moving_mat), diag(moving_mat_simul))
    inflows = zeros(50)
    inflows_simul = zeros(50)
    for i = 1:50
        inflows[i] = (sum(moving_mat[:,i])-diag(moving_mat)[i])/50
        inflows_simul[i] = (sum(moving_mat_simul[:,i])-diag(moving_mat_simul)[i])/50
    end
    mig_corr_in = cor(inflows, inflows_simul)

    err_vec[64] = (1-mig_corr)
    err_vec[65] = (1-mig_corr_in)
    =#

    #####################Additional Collge Moments################
    coll_moments_additional = [0.344, 0.162, 0.265, 0.4017, 0.5476, 0.0689, 0.1820, 0.300, 0.473, 0.648, 0.258, 0.636]
    coll_overall = mean(data_simul[:,13], weights(weights_all))
    println("OVERALL COLLEGE ATTENDANCE: ", coll_overall)
    err_vec[63] = (coll_overall - coll_moments_additional[1])/coll_weight
    index_moment_coll = 2
    index_err_vec = 64

    #by income quartile
    for i = 1:4
        minpctile = (i-1)*25
        maxpctile = (i)*25
        restrict_pctile_pt1 = data_simul[data_simul[:,ncol-1].>=minpctile, :]
        restrict_pctile_pt2 = restrict_pctile_pt1[restrict_pctile_pt1[:,ncol-1].<maxpctile, :]
        weights_temp = fweights(restrict_pctile_pt2[:,14])
        coll_pctile = mean(restrict_pctile_pt2[:,13], weights(weights_temp))
        err_vec[index_err_vec] = ((coll_pctile - coll_moments_additional[index_moment_coll])/coll_weight) / 2
        index_err_vec+=1
        index_moment_coll += 1
        println("College Attendance, Quartile ", i, ": ", coll_pctile)
    end

    #by ability
    for i = 1:5
        restrict_ability = data_simul[data_simul[:,11].==i, :]
        weights_temp = fweights(restrict_ability[:,14])
        coll_abil = mean(restrict_ability[:,13], weights(weights_temp))
        err_vec[index_err_vec] = ((coll_abil - coll_moments_additional[index_moment_coll])/coll_weight) / sqrt(5)
        index_err_vec+=1
        index_moment_coll += 1
        println("College Attendance, Ability Level ", i, ": ", coll_abil)
    end

    #by parent education attainment status
    for i = 1:2
        restrict_pcoll = data_simul[data_simul[:,12].==i,:]
        weights_temp = fweights(restrict_pcoll[:,14])
        coll_pcoll = mean(restrict_pcoll[:,13], weights(weights_temp))
        err_vec[index_err_vec] = ((coll_pcoll - coll_moments_additional[index_moment_coll])/coll_weight) / sqrt(2)
        index_err_vec+=1
        index_moment_coll += 1
        println("College Attendance, Parent Coll: ", i, ":", coll_pcoll)
    end


    index_err_vec = 75

    #75-124: outflow
    for l = 1:prim_2000.nl
        mig_out, mig_out_simul = diag(moving_mat)[l], diag(moving_mat_simul)[l]
        err_vec[index_err_vec] = ((mig_out - mig_out_simul)/mig_out)/sqrt(prim_2000.nl)
        index_err_vec += 1
    end

    #125-174: inflow (where things might get a bit ugly)
    state_inflow = package_moments[5]
    restrict_move = data_simul[data_simul[:,1].!=data_simul[:,2],:] #movers in simulated data
    n_move = sum(restrict_move[:,14]) #total weight of movers
    for l = 1:prim_2000.nl
        mig_in = state_inflow[l]
        restrict_move_state = restrict_move[restrict_move[:,2].==prim_2000.fips[l], :]
        n_outflow = sum(restrict_move_state[:,14])
        mig_in_simul = n_outflow/n_move
        err_vec[index_err_vec] = ((mig_in - mig_in_simul)/mig_in)/sqrt(prim_2000.nl)
        index_err_vec += 1
    end

    #175-224: iim
    state_iim = package_moments[6].*100
    ncols = length(data_simul[1,:]) #useful for later
    restrict_median = data_simul[data_simul[:,ncols-1].<50, :]

    #construct state means
    for l = 1:prim_2000.nl
        state_iim_simul = mean(restrict_median[:,ncols][restrict_median[:,1].==prim_2000.fips[l]])
        err_vec[index_err_vec] = ((state_iim[l] - state_iim_simul)/state_iim[l])/sqrt(prim_2000.nl)
        index_err_vec += 1
    end

    ######race moments
    index_err_vec = 225

    #get racial weights
    data_simul_w = data_simul[data_simul[:,20].==1, :]
    data_simul_b = data_simul[data_simul[:,20].==2, :]
    data_simul_h = data_simul[data_simul[:,20].==3, :]
    weights_b = fweights(data_simul_b[:,14])
    weights_h = fweights(data_simul_h[:,14])
    weights_w = fweights(data_simul_w[:,14])

    #college attendance
    coll_simul_b = mean(data_simul_b[:,13], weights(weights_b))
    coll_simul_h = mean(data_simul_h[:,13], weights(weights_h))
    err_vec[225] = ((coll_simul_b - moments_race[1])/moments_race[1])
    err_vec[226] = ((coll_simul_h - moments_race[2])/moments_race[2])

    #college attendance by region
    data_simul_r_1 = data_simul[data_simul[:,21].==1, :]
    data_simul_r_2 = data_simul[data_simul[:,21].==2, :]
    data_simul_r_3 = data_simul[data_simul[:,21].==3, :]
    data_simul_r_4 = data_simul[data_simul[:,21].==4, :]
    weights_r_1 = fweights(data_simul_r_1[:, 14])
    weights_r_2 = fweights(data_simul_r_2[:, 14])
    weights_r_3 = fweights(data_simul_r_3[:, 14])
    weights_r_4 = fweights(data_simul_r_4[:, 14])

    coll_r_1 = mean(data_simul_r_1[:, 13], weights(weights_r_1))
    coll_r_2 = mean(data_simul_r_2[:, 13], weights(weights_r_2))
    coll_r_3 = mean(data_simul_r_3[:, 13], weights(weights_r_3))
    coll_r_4 = mean(data_simul_r_4[:, 13], weights(weights_r_4))
    err_vec[227] = ((coll_r_1 - moments_race[3])/moments_race[3])
    err_vec[228] = ((coll_r_2 - moments_race[4])/moments_race[4])
    err_vec[229] = ((coll_r_3 - moments_race[5])/moments_race[5])
    err_vec[230] = ((coll_r_4 - moments_race[6])/moments_race[6])

    #migration
    mig_simul_b = mean(data_simul_b[:,15], weights(weights_b))
    mig_simul_h = mean(data_simul_h[:,15], weights(weights_h))
    err_vec[231] = ((mig_simul_b - moments_race[7])/moments_race[7])
    err_vec[232] = ((mig_simul_h - moments_race[8])/moments_race[8])

    #period 3 earnings
    earn_simul_b = mean(data_simul_b[:,9], weights(weights_b))
    earn_simul_h = mean(data_simul_h[:,9], weights(weights_h))
    earn_simul_w = mean(data_simul_w[:,9], weights(weights_w))
    err_vec[233] = (((earn_simul_b/earn_simul_w) - moments_race[9])/moments_race[9])
    err_vec[234] = (((earn_simul_h/earn_simul_w) - moments_race[10])/moments_race[10])

    #######optional spitting out of values#######
    if report == 1

        println("Rank-rank ige: ", round(ige, digits = 3))
        println("Mean Young Earnings, HS: ", round(mean_earn_hs_p2, digits = 3))
        println("STD Young Earnings, HS: ", round(std_earn_hs_p2, digits = 3))
        println("Mean Earnings Growth, HS: ", round(mean_earn_hs_growth, digits = 3))
        println("STD Earnings Growth, HS: ", round(std_earn_hs_growth, digits = 3))
        println("Mean Old Earnings, HS: ", round(mean_earn_hs_p3, digits = 3))
        println("STD Old Earnings, HS: ", round(std_earn_hs_p3, digits = 3))

        println("Mean Young Earnings, COLL: ", round(mean_earn_coll_p2, digits = 3))
        println("STD Young Earnings, COLL: ", round(std_earn_coll_p2, digits = 3))
        println("Mean Earnings Growth, COLL: ", round(mean_earn_coll_growth, digits = 3))
        println("STD Earnings Growth, COLL: ", round(std_earn_coll_growth, digits = 3))
        println("Mean Old Earnings, COLL: ", round(mean_earn_coll_p3, digits = 3))
        println("STD Old Earnings, COLL: ", round(std_earn_coll_p3, digits = 3))

        println("HS Mig Rate: ", round(mig_hs, digits=3))
        println("COLL Mig Rate: ", round(mig_coll, digits=3))

        println("Move-stayer HC difference, HS: ", round(hc_move_hs-hc_stay_hs, digits = 3))
        println("Move-stayer HC difference, COLL: ", round(hc_move_coll-hc_stay_coll, digits = 3))


        #println("OVERALL MIGRATION CORRELATION: ", mig_corr)
        #println("OVERALL MIGRATION INFLOW CORRELATION: ", mig_corr_in)
        println("Average Staying Rate: ", round(mig_avg, digits = 3))
        println("Staying Rate STD: ", round(mig_stdev, digits = 3))
        println("Population-Moving Correlation: ", round(corr_pop_mig_in, digits = 3))
        println("Population-Staying Correlation: ", round(corr_pop_stay, digits = 3))
        println("Share of moves to Proximal Locations: ", round(prox_share, digits = 3))
        println("Period-3 Migration Rate: ", round(mig_p2_rate, digits=3))
        #println("Period-3 Home Migration Rate: ", round(mig_home_rate, digits=3))

        println("Average Time Investment: ", round(time_inv_avg, digits = 3))

        #racial moments
        println("College Attendance, Black: ", round(coll_simul_b, digits=3))
        println("College Attendance, Hispan: ", round(coll_simul_h, digits=3))
        println("College Attendance, Region 1: ", round(coll_r_1, digits=3))
        println("College Attendance, Region 2: ", round(coll_r_2, digits=3))
        println("College Attendance, Region 3: ", round(coll_r_3, digits=3))
        println("College Attendance, Region 4: ", round(coll_r_4, digits=3))
        println("Migration, Black: ", round(mig_simul_b, digits=3))
        println("Migration, Hispan: ", round(mig_simul_h, digits=3))
        println("Earnings, Black: ", round(earn_simul_b/earn_simul_w, digits=3))
        println("Earnings, Hispan: ", round(earn_simul_h/earn_simul_w, digits=3))
    end

    err_vec #return estimation error
end

#function to optimize in estimation of model
function Objective_function(guess::Array{Float64,1}, package::Vector{Matrix{Any}}, package_moments::Vector{Matrix{Float64}};
    nsims::Int64 = 2000, ret::Int64 = 0, report::Int64 = 1, seed::Int64 = 1234, cfact::Int64 = 0, l_select::Int64 = 0, aopt::Int64=0, statcheck::Int64 = 0, short::Int64=0)

    #report current guess
    println("\n")
    println("CURRENT GUESS")
    println(round.(guess, digits = 3))

    #preallocation
    data_simul = zeros(2,2)
    data_simul_p25 = zeros(2,2)
    output_simul = [0, 0, 0, 0, 0]
    err_vec = zeros(234)
    error = 0.0
    moments_targ = package_moments[1]
    flag_stop = 0

    #check for flipping kill switch
    if guess[1]>1.0 || guess[1]<0.0 || guess[2]>=1.0 || guess[2]<=-1.0 || guess[4]<=0.0
        flag_stop=1
    elseif guess[5]<=0 || guess[6]<=0.0 || guess[6]>=1.0 || guess[7]<=0.0 || guess[7]>=1.0 || guess[8]<=0.0
        flag_stop=1
    elseif guess[14]<=0.0 || guess[18]<=0.0
        flag_stop=1
    end

    #kill if bad paramter values
    if flag_stop == 1
        error = Inf
    else

        prim_2000, res_2000 = Initialize(package, package_moments)
        prim_2010, res_2010 = Initialize(package, package_moments; p2 = 1)

        println("Solving Model for 2000 . . .")
        grids, est = Initialize_grids(prim_2000, guess) #initialize computing grids
        V_iterate(prim_2000, est, res_2000, grids; l_select = l_select, cfact=cfact, aopt=aopt) #compute policy and value function.

        println("Solving Model for 2010 . . .")
        V_iterate(prim_2010, est, res_2010, grids; l_select = l_select, cfact=cfact, aopt=aopt) #compute policy and value function.

        #simulations
        println("Simulating Data . . .")
        output_simul = Simulate(prim_2000, prim_2010, est, res_2000, res_2010, grids, package; nsims = nsims, seed=seed, aopt=aopt, statcheck=statcheck, short=short) #simulate CHKS cohort
        err_vec = Compute_error(prim_2000, output_simul, package_moments; report = report)

        #report and adjust error
        W = package[14]
        error = err_vec' * W * err_vec

        ###opti;onal reporting
        #if report == 1
            #println("Marriage-IIM State-Level Correlation: ", iim_corr)
            #println("IIM Mean: ", mean(state_iim_simul[:,1]))
            #println("IIM SD: ", std(state_iim_simul[:,1]))
        #end

        #in estimation mode
        if ret == 0
            finalize(res_2000)
            finalize(res_2010)
            finalize(output_simul)
            @everywhere GC.gc()
        end
    end

    println("Current error: ", error)
    println("\n")


    #write to log file
    open("logfile.txt", "a") do f
        print(f, round.(guess, digits = 3))
        print(f, "\n")
        print(f, error)
        print(f, "\n")
        print(f, "\n")
    end


    #return deliverables
    deliverables = error

    if ret == 1
        deliverables = [error, output_simul[1], data_simul_p25, output_simul[4], err_vec]
    end

    deliverables
end
